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integral decomposition and a worm operator which is local in imaginary time. It generates states with a fixed 
number of particles and respects other exact symmetries. Observables like the equal-time Green's function can 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) simulation is a powerful and 
versatile method for the investigation of thermodynamic prop- 
erties of many-body systems. When generating a Markov- 
chain of configurations using a Metropolis scheme |1J], it is 
known that updates based on local changes are inefficient, par- 
ticularly near critical points. At transitional points this type 
of algorithm gives very large autocorrelation times, a phe- 
nomenon one refers to as 'critical slowing down' Q|. By 
developing non-local update schemes, this problem has been 
overcome for second order phase transitions. The Wolff single 
cluster algorithm |3] and the Swendsen-Wang multiple clus- 
ter algorithm |4], both used to solve classical physics prob- 
lems, were successful applications of this idea. In the same 
spirit, loop algorithms |5j, |a |7J] were developed for the study 
of discrete quantum systems. New configurations are obtained 
by flipping clusters in the form of loops. The systematic er- 
ror caused by the the Suzuki-Trotter approximation was 
eliminated by formulating the world-line algorithms directly 
in continuous imaginary time lUol Hill . In the worm algo- 
rithm 1 12], the partition function is embedded in an extended 
configuration space, allowing a direct and exact evaluation 
of the one-body Green's function. The concept of non-local 
loop updates has also been implemented in the Stochastic Se- 
ries Expansion (SSE) representation lfl3ll . leading to 'opera- 
tor loop' update algorithms [14, El and 'directed loop' al- 
gorithms llfjl 1 1711 . which are a further optimization of the 
loop construction. The general idea behind this is to construct 
moves in a locally optimal way fiUl . 

All the non-local update world-line algorithms which are 
mentioned above sample the grand-canonical ensemble. In 
this way one can generate configurations with e.g. varying 
magnetization or occupation number. Results for the canon- 
ical ensemble are obtained by using only the configurations 
with the right particle number [19] or by rejecting loop up- 
dates which change this number |6]. It is clear that this is an 
inefficient way of working. Sampling the canonical ensemble 
directly would be advantageous when studying systems where 
particle-conservation is important. One example is the transi- 
tion between the superfluid and the Mott phase in the Bose- 
Hubbard model at constant filling. This transition belongs to 



the (d+l)-dimensional XY universality class 1 20] . When en- 
tering the superfluid phase, it becomes difficult to keep the 
number of bosons constant and tuning the chemical potential 
can become a very time-consuming task. When simulating 
mesoscopic systems like superconducting grains 12 U 12211 or 
atomic nuclei l23l l24ll . it is primordial to keep the number 
of particles constant. A world-line algorithm satisfying these 
conditions is presented 12511 and discussed in detail in this pa- 
per. Besides particle-number conservation, the algorithm al- 
lows to include other symmetry-projections. It is constructed 
in such a way that all moves are accepted, which makes it effi- 
cient to run and easier to code. Though working in the canoni- 
cal ensemble, our algorithm is still able to generate configura- 
tions with different winding numbers, in contrast to the local 
world-line method by Hirsch et al. [26]. We will demonstrate 
the versatility of the method by applying it to bosons and to 
paired fermions. 

II. THE ALGORITHM 

Practically all QMC methods are based on a decomposition 
of the evolution operator e~$ H . The trick is to break up this 
operator into pieces which can be handled exactly 1 27] . Gen- 
erally one can write the Hamiltonian as consisting of an easy 
part Ho and a residual interaction V, 

H = H -V. (1) 

The minus sign has been included to ease notations further on. 
For such a Hamiltonian, one can make an exact perturbative 
expansion in V of the evolution operator, using the following 
integral expression Hill : 

e-^ H = V f dt m f m dt m ^--- f 2 dt\V(h) 
„^ Jo Jo Jo 

V(t2)---V{t m )e-^, (2) 

with V(t) =exp(—tHo)Vexp(tHo) and (3 the inverse temper- 
ature (also called imaginary time). The basic idea of the 
continuous-time loop algorithm |6, 10] and the worm algo- 
rithm 1 12] is to insert two adjoint world-line discontinuities. 
By propagating one of these discontinuities (which are called 
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the mobile 'worm head' and stationary 'worm tail') through 
lattice space and imaginary time, the configuration changes in 
such a way that detailed balance is fulfilled. At that point one 
is not sampling according to the partition function Tr(e~^ H ), 
but according to an extension hereof, 

Tr(W^e- zH We-^ H ), (3) 

with x the imaginary time interval between the worm 'head' 
operator and 'tail' operator W. The worm head can be 
creating or annihilating, depending on the choice of W. Af- 
ter some Markov steps, the worm head bites its tail and the 
discontinuities are removed. Only configurations with con- 
tinuous world-lines can contribute to the statistics according 
to Tr(e~^ H ). In contrast to these algorithms, we will work 
with a worm which is local in imaginary time. The evolution 
operator extended by such a local worm (an imaginary time- 
independent operator A to be defined below) reads 

U / ($,l)=e- xH Ae-^ H , (4) 

where x can be regarded as the imaginary time at which the 
worm operator is inserted. We will show that by working 
with a local worm operator one can construct a very efficient 
sampling method, provided that the worm operator commutes 
with the residual interaction (AV = VA). If A furthermore 
commutes with the generators of a symmetry of Ho and V, 
one can restrict the sampling to configurations with those spe- 
cific symmetries, leading to symmetry-projected results. In 
particular one can sample the canonical ensemble with a worm 
operator that conserves particle number. We would like to em- 
phasize at this point that the algorithm stated below does not 
depend on the specific structure of A. The operator A has to 
be chosen in such a way that an ergodic Markov chain can be 
constructed, and therefore it will depend on the specific form 
of the interaction V. 

If one decomposes the trace (restricted to the wanted parti- 
cle number and symmetry) of £/'(p,x) using Eq. (|2ji and in- 
serts complete sets of eigenstates of Ho at all imaginary times, 
one obtains a set of integrals which can be evaluated using 
Monte Carlo sampling. The Markov process will sample the 
configurations proportional to the weights 

W(m,i,t,x) = {h\V\h)e-^-' l)E n (i 1 \V |/ 2 }^ (,3 ^ 2)£ '2 • ■ • 
(k-i \V\i L )e-^ E 'L (i L \A\i R )e-^ E 'R (i R \V\i R+l ) ■ ■ ■ 
(imr-l \V\i m )e-^-^ E >- (i m \V\i Q )e-^-^ E -o. (5) 

Each configuration is specified by an order m (the num- 
ber of interactions), a set i of eigenstates of Ho (with 
i a shorthand notation for all the intermediate states 
|z'o), . . . , \i R ), | /,«)), interaction times t\, . . . ,t m , and the 
worm insertion time x. We use the notation E{, = (ij\Ho\ij). 
The configuration to the left of the worm is changed by 
the worm operator into the configuration | i R ). We use the sub- 
script L (R = L+ 1) to indicate the eigenstate ( \i R )) and 
interaction time ti (t R ) just before (after) the worm operator 
in imaginary time. We will call the configurations for which 
il = i R diagonal configurations. By choosing the worm opera- 
tor such that its diagonal elements are constant (i.e. (i|A|i) = c 



for all basis states |/)), the sum of the weights of all diago- 
nal configurations is proportional to the particle-number pro- 
jected trace of the evolution operator t/(P). This is nothing 
else than the canonical partition function Ttff(e~^ H ), with Tr^ 
the particle-number projected trace. Hence, sampling the con- 
figurations proportional to their weights W(m, i,t,x) leads to 
a sampling of the canonical ensemble. The Markov process is 
set up using the Metropolis-Hastings algorithm 1 1, 28], hereby 
sampling in an extended space according to Trjv[I/'(P,x)]. At 
each Markov step only a few of the factors of Eq. (|5} are al- 
tered by the worm operator which moves to a new point in 
imaginary time. These worm operator moves will be con- 
structed in a such a way that detailed balance is fulfilled 
locally at each Markov step. Therefore it is also fulfilled 
when going from one diagonal configuration to another. It 
takes a number of Markov steps before diagonal observables 
(i.e. observables which commute with Ho) can be measured 
again. While each Markov step contains only local changes, 
the chain of steps between two diagonal configurations cor- 
responds to a global update. Non-diagonal operators can be 
measured using the fact that one samples in an extended space. 
By keeping track of the worm moves between two diagonal 
configurations, statistics for the expectation values of non- 
diagonal operators can be collected, similar to the way one 
evaluates the one-body Green's function in the worm algo- 
rithm 111211 . Our method however will lead to much better 
statistics for equal-time non-diagonal operators, because the 
worm operator is always local in imaginary time. 

Before shifting the worm operator over some imaginary 
time interval Ax, a direction D has to be chosen. One can 
choose between the directions D = R (higher values of x) and 
D = L (lower values of x) with some probability P(D), to be 
specified later. The presence of the exponentials in Eq. Q 
inspires us to choose the time shift Ax proportional to an ex- 
ponential distribution, 

P(Ax)dAz = e D e- £DAT dAz, (6) 

with to an optimization parameter. In shifting the worm op- 
erator from x to a new imaginary time x' = x + Ax, the worm 
operator can encounter an interaction operator V at some time 
in between. Assume the direction R is chosen and the worm 
operator meets an interaction at time t R . Let us first consider 
the situation where the worm operator moves through this in- 
teraction, without annihilating it, 

<iL|A|i*)<i*|V|i*+i) — > (i L \V\i' R ){i' R \A\i R+l ). (7) 

When passing the interaction, the intermediate state can 
change. A convenient way to pick one of these possible 
changes is to choose the new configuration proportional to its 
weight 

(/ L |y|4)(4|A|/fl +1 ) 

r R +i,L{i R ) = j. ; ,. r ■ (») 

(i L \VA\i R+ i) 

Part (a) of Figure ^ shows a diagrammatic representation of 
the different ways in which a general one-body worm operator 
A = jCijaJaj (for some constants cy ) can pass a similar 
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interaction V. The worm operator is represented by a curly 
line and the interaction by a vertical straight line. For this type 
of worm and interaction operator there are always at most four 
ways in which the intermediate state can change. It should be 
noted that our choice Eq. (|8} is not unique and possibly more 
optimal choices can be found [18]. Because of this choice 
however, there appears a factor 



of interactions, 



(iL\V\i' R )(i' R \A\i R+l )P L , R +i(iR) 
(i L \A \i R )(i R \V\ i R + i)Pr+\ , L (i R ) 
{i L \VA\i R+l } 
(fe|AV|iji + i) ' 



(9) 



in the acceptance factor of the Metropolis-Hastings algorithm. 
Every time the worm operator passes an interaction an anal- 
ogous factor appears, depending on the direction D of propa- 
gation. Therefore it's advantageous to impose on A the condi- 
tion, 



AV-VA = 0, 



(10) 



because then n DD i = 1 in all cases, and we do not have to 
worry about this normalization factor anymore. Furthermore, 
Eq. J ensures that the worm operator can always pass the 
interaction it encounters. If one would choose a worm opera- 
tor A that does not satisfy this condition, as in grand-canonical 
algorithms, it is possible the worm operator cannot pass the 
interaction and the direction of propagation has to change, in 
that way undoing changes previously made. It is intuitively 
clear that these 'bounc es' give rise to a slow decorrelation and 
should be avoided Ilafl7lfl8ll . In the directed loop algorithm 
one increases the efficiency of the loop update by minimiz- 
ing the appearance of bounces, but they cannot be eliminated 
completely because of the reversibility condition. In what fol- 
lows we will assume the condition Eq. dlOt is fulfilled, mak- 
ing the algorithm bounce-free. We will drop the factors n DD i 
to ease the equations. After passing through the interaction 
at time to, one has to choose a new imaginary-time shift At. 
However, one can avoid generating a new random number by 
taking the new shift as follows: 



At = (Y - t D ) 



old 



(ID 



where the parameter to has been updated after passing the 
interaction. 

The choice of the parameters £o follows from detailed bal- 
ance. Because the time shifts At of the worm operator are 
chosen by Eqs. and il 1> . adding the constraint 



E R -E L = s L -£ R , 



(12) 



ensures that all the exponentials which appear in the ac- 
ceptance factor of the Metropolis-Hastings algorithm cancel, 
leading to an efficient sampling method. So in practice one 
can choose any positive value for £l and e R , as long as Eq. 
d!2i is fulfilled at each step of the worm movement. To con- 
clude, we write down the acceptance factor for the above pro- 
cedure when the worm operator does not change the number 



W(m,i',t',x')P(i',t',x' ^i,t,x) 
W(m,i,t,z)P(i,t,x— y i',f',x') 

(£/)') initial 



(13) 



final 



where (Ez))nnai ((ED')initiai) is the value of Eo (£/)') at the en d 
(beginning) of the worm operator move into direction D, and 
D' denotes the opposite direction of D. The actual acceptance 
probability is given by min(l , q), according to the Metropolis- 
Hastings algorithm. 

Let us now introduce a number of steps, which allow to 
change the number of interactions in a reversible way. We 
want the acceptance factor of such updates to be local, i.e. the 
probability to pass, create or annihilate an interaction should 
only depend on the properties of the state at that point in 
imaginary space-time. We consider three extensions of the 
procedure outlined above where no interactions are created or 
deleted: 

• At the beginning of the Markov step, we introduce the 
possibility that the worm operator creates a new inter- 
action with probability co, which depends on the direc- 
tion D of propagation. This creation will also change 
the intermediate state: 



(k\A\i R ) — (iL\V\i')(i'\A\i R ) 



(14) 



assuming again the worm operator is moving in the 
D = R direction. The new intermediate state |/') will 
be chosen with probabilities 



PRLif) 



(iL\V\i')(i'\A\i R ) 
(i L \VA\i R ) 



(15) 



For a worm operator move in the D = L direction, prob- 
abilities Pw{i') can be defined in an analogous way. 
Figure^(parts (b) and (c)) shows a diagrammatic rep- 
resentation of the insertion of a one-body interaction at 
the beginning of the worm move. For a diagonal con- 
figuration only the diagonal part of A contributes to the 
matrix element (ii\A\i R ). In this situation the worm op- 
erator is represented by little circles and all world-lines 
are continuous. We will call this the 'diagonal worm'. 

• When the worm operator arrives at an interaction, one 
also has to consider the possibility of annihilating that 
interaction. Suppose the interaction can be deleted. Let 
ao be the probability to annihilate the interaction while 
continuing the worm update, and so the probability to 
annihilate the interaction and stop the worm update. 
Then 1 — ao — so is the probability to pass through that 
interaction and continue the worm operator. 

• To maintain reversibility, one also has to include the 
possibility that the construction of the Markov step does 
not halt at the moment the worm operator has finished 
a shift At without encountering an interaction. At that 
point one has to choose between stopping the worm op- 
erator, or to continue, with the possibility of inserting a 
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new interaction at that point. Let fo be the probability 
to continue the worm operator without inserting an in- 
teraction, go the probability to insert an interaction and 
to continue the worm operator, then I — fo — gD will 
be the probability to stop the worm operator, without 
inserting an interaction. 

After creating, annihilating or passing an interaction, a new 
time shift Ax should again be chosen according to Eqs. (|6j 
or (1111 . Note that the parameter fo is redundant: jumping 
with a parameter Ed and continuing the worm operator unal- 
tered with probability fo is statistically equivalent to making 
ajump with parameter £d(1 —/d), and then choosing between 
either stopping the worm operator or inserting an interaction 
and move on. Therefore, one can set /d — without loss 
of generality. We now determine how the other parameters 
should be chosen in order to satisfy detailed balance. Assume 
a direction D is chosen. When no interaction is inserted at the 
beginning of the worm move, a factor 



where 



o _ £ £>'(1 -gzy. 

1 - CD 



(16) 



appears in the acceptance factor. If on the other hand an inter- 
action is created at the initial time x of the worm operator, this 
will lead to a factor 



with 



Id 



Hdd< 



HDD'S pi 
CD 



{i D \VA\ijy) 



(17) 



(18) 



{h\A\i D ') ' 

A new intermediate state is chosen with probabilities Eq. Jl 5I >. 
At the end of the Markov step, the worm operator will annihi- 
late an interaction or not, leading to extra factors in the global 
acceptance factor which have the inverse form of Eqs. (1171 
and d!6l >. because of the inverse symmetry between beginning 
and end of the move. At intermediate points, we can encounter 
the following situations. The worm operator can stop after a 
shift Ax between two interactions, insert an interaction and 
move on. The inverse situation of this can also occur, when 
an interaction is annihilated and the worm operator moves on. 
Or the worm operator can simply pass an interaction without 
annihilating it. In order to have a good total acceptance fac- 
tor, we will require that these intermediate steps do not con- 
tribute to the acceptance factor. This condition leads to the 
constraints 



HdD'Od' = £DgD, 

aD + SD = a D i+s D i- 



(19) 
(20) 



Apart from that, we want the sampling to be as uniform as 
possible, which suggests the condition q D — q c D . Putting all 
this together, the total acceptance factor is given by 



W(rri ,i' ,f ,?)P(tri ,i' ,t' ,x' -> m,i,t,x) 
W(m,i,t,x)P(m,i,t,x — > m',i',t',x') 

(go) initial 
(<?£>') final ' 



(21) 



qD = £d> + *Cdd> ( s d' -a D )- 



(22) 



The factor (gD)mitiai/finai has to be evaluated at the beginning 
and the end of the Markov step in direction D (with D' the 
opposite of D). The creation probability is given by Eq. Mil . 



CD 



HDD'S pi 

qD 



(23) 



We still have to determine how to choose the direction D. The 
acceptance ratio of Eq. Mil inspires us to choose the direction 
of the move with probabilities 



P(D = R) = 



qR 
Rlr ' 
P(D = L) = 

Rlr 



(24) 
(25) 



with Rlr = qR + qL- By accepting all worm operator moves a 
distribution given by 



W'{m,i,t,z) =RLRW(m,i,t,z), 



(26) 



will be sampled, instead of the distribution W(m,i,t,x). Be- 
cause the factors Rlr fluctuate only mildly in practice, accept- 
ing all moves still leads to a a very useful sampling method. 
It speeds up the algorithm and reduces the complexity of the 
code. Finite temperature observables can still be evaluated by 
taking the extra weighting factor into account: 



<e) P 



Tr[e 



-P(*o-v; 



Q] 



7>[ e -P(flb-v)] 



(27) 



Each time the worm operator creates, annihilates or passes 
an interaction, the parameters Ed, cd, cid, sd, go are deter- 
mined by Eqs. (I12> . &20\ , i22i and ( 1231 . This still leaves a lot 
of freedom. We will focus here on two limiting cases. First 
we will consider the case where one of the two parameters 
£r and El is always zero. In that way it can occur the time 
shift Ax becomes infinite. This amounts to the worm opera- 
tor directly jumping to the next interaction, which speeds up 
decorrelation in imaginary time direction. In order to obtain a 
worm move that changes the configurations as much as possi- 
ble, the parameters gL, gR, «l and cir are maximized. The set 
of parameters obtained in this way is shown in Table|l]for the 
case Er > El- We will call this solution A. The solution for 
the case El > Er can be found by interchanging the subscripts 
L and R. Note that in this solution the worm operator always 
starts to move into the direction of the highest diagonal en- 
ergy. It is clear that whenever the worm operator is moving 
in the direction of the highest diagonal energy or whenever 
Er = El, the time shift Ax becomes infinite. There are a num- 
ber of extra conditions one should keep in mind. Assume the 
worm operator starts to move in the direction D = R (because 
Er > El). When the worm operator arrives at an interaction 
that can be annihilated, one has to determine El, Er and 9\[lr 
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FIG. 1 : A diagrammatic representation for one-body worm operator 
moves. The worm operator is represented by a curly line and the 
interaction V by a solid vertical line. We distinguish between the 
following updates, (a) When the worm operator moves in the D = R 
direction, it can encounter some interaction. The worm operator can 
pass the interaction, in that way changing the intermediate state. For 
general one-body worm and interaction operators, there are at most 
four possible ways in doing this, (b)-(c) At the beginning of the 
worm move, we introduce the possibility of inserting an interaction. 
When the initial worm is diagonal (represented by circles), a number 
of interaction insertions of the type shown in (b) are possible. In 
part (c) the initial worm operator is not diagonal, and an interaction 
insertion can make the worm diagonal or not. 



after the annihilation. If Er > El is still satisfied then, sr and 
Or from Table U are the correct probabilities to stop or con- 
tinue the worm operator. If now El > Er on the other hand, 
one has to use the solution sr = min(l, ]C R ) and cir = 0, but 

the worm operator keeps moving in the same direction. The 
time shift of the worm operator is only finite when it moves in 
the direction D and Ed < E' D . Note also that only gL is men- 
tioned in Table|l] because the parameter go has only meaning 
when the time shift is finite. In the present solution however, 
a problem arises whenever El = Er. In this case Er = El = 
and the time shift At is always infinite. Because sr = sl = 
in addition, the worm operator never halts. As a consequence 
configurations with a diagonal worm will never be sampled. 
This can be solved by proposing a small but finite stopping 
probability. This alternative solution for the case El = Er is 
also given in Table|l] The global parameter should be taken 
small (such that < < "JioD' f° r a U diagonal configurations) 
but not zero, and can be chosen in order to optimize the decor- 
relation between successive evaluations of observables. Note 
Rlr of Eq. i26\ takes the constant value 20. 

Another possibility to find algorithm parameters follows 
from the idea that we want the step size At to be of the order 
of the mean imaginary time interval between two interactions. 
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TABLE I: A set of algorithm parameters satisfying Eqs. <12> . 1201 . 
1221 and 1231 for the cases El = Er and El < Er (otherwise inter- 
change the subscripts L and R). We call this solution A, for which 
one of the parameters is always zero. 
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TABLE II: An alternative set of parameters for which one of the 
parameters is always 1{lr- We call this solution B. The parameter 
can be chosen to optimize the algorithm. 



Therefore we consider the case where one of the two param- 
eters Er and El is 9{lr- As a consequence the time shift is 
always finite. For Er > El such a set of parameters is given in 
Table |ll| Again, the case El = Er needs an alternative solu- 
tion, since otherwise Rlr would be zero for diagonal configu- 
rations. We will refer to this solution as solution B. 

In short, the algorithm is based on a time-dependent per- 
turbation in the interaction V (see Eq. (0), so there is no 
systematic error arising from time discretization. Because we 
choose time shifts of a worm operator according to Eq. (|6jl 
the diagonal part of the Hamiltonian is handled exactly. There 
are a number of algorithm parameters, which have to satisfy 
Eqs. d!2l >. J20t . i22i and ( I23l l. We have derived two sets of 
algorithm parameters, satisfying these equations. In the first 
set (solution A) one of the parameters Ed is always zero and 
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in the second (solution B) it is equal to 9{lr- Therefore the 
main difference between these two solutions will be the size 
of the imaginary time shift Ax. Other algorithms where 8« 
or El take values between and s\£l« can be constructed in a 
similar way, taking for now only these limiting cases. In the 
next section our QMC algorithm will be applied to the Bose- 
Hubbard model. The efficiency of the two different solutions 
leading to different algorithms will be compared in this con- 
text. 



III. APPLICATION TO THE BOSE-HUBBARD MODEL 

Ultracold bosonic atoms in an optical lattice are described 
by the Bose-Hubbard model I29ll30ll3lll . 



H = -t£b f i bj + -'£n i (n i -l), 



(28) 



with b] (bj) the boson creation (annihilation) operator on site 
i, tii the number operator on site i and denoting nearest 
neighbors. The lattice has N s sites, occupied by N bosons. 
The parameter t is the tunneling amplitude and U is the on- 
site repulsion strength. We will restrict the discussion here 
to the one-dimensional homogeneous case without trap. At 
low values of U/t the system forms a compressible super- 
fluid. This phase is characterized by a gapless excitation spec- 
trum and long-range phase coherence. By increasing U/t, a 
quantum phase transition from a superfluid state to a Mott 
insulating state is achieved at integer densities. In the pure 
Mott phase the bosons are localized at individual lattice sites 
and all phase coherence is lost due to quantum fluctuations. 
In addition, density fluctuations disappear when entering the 
Mott phase and a gap appears in the excitation spectrum. This 
phase driven transition belongs to the Berezinskii-Kosterlitz- 
Thouless (BKT) 13 2l 13311 universality class in one dimension. 
The Bose-Hubbard Hamiltonian can be rewritten in the form 
Eq. Q, 

U 



Ho = — ^«i(n ; - 1), 

L i 

V = t^bjbj. 



(29) 



As argued above, it is advantageous to take the worm operator 
A such that it commutes with V. An operator that satisfies this 
condition is given by 



A = — Y Hi - 



(30) 



with N a c-number to be optimized. In our calculations this 
parameter is always set equal to the total number of particles. 
We have checked our code by comparing with exact diagonal- 
ization results for small lattices. Ergodicity was tested numer- 
ically. Power law behavior of correlation functions coincides 
with predictions from bosonization theory for large lattices. 



The one-body Green's function G(r) — (b\bi +r ) is calcu- 
lated by updating the entry r in a histogram for G(r) at ev- 
ery Markov step. The function G(r) can be normalized di- 
rectly from the diagonal / non-diagonal worm fraction. The 
non-diagonal worm components b]bi+ r of Eq. ( I30i can be 
given a different weight, leading to a worm matrix represen- 
tation of the symmetric Toeplitz type (i.e. a symmetric matrix 
with constant values along negative-sloping diagonals). Such 
a worm operator still commutes with the interaction part V of 
the Hamiltonian. By giving some worm components a bigger 
weight, the corresponding components G(r) will be updated 
more often, leading to a higher accuracy and mimicking flat 
histograms 1 34 ] . The condensed fraction p 6 can be determined 
via 



(31) 



The superfluid fraction can be determined using the winding 
number 1 3f 



(W 2 )N 2 
2tN$ 



(32) 



where (W 2 ) is the mean square of the winding number oper- 
ator in one dimension. Figure [2] shows the condensed and su- 
perfluid fraction for a uniform one-dimensional system of 128 
sites at a density of exactly one particle per site. Calculations 
were performed at an inverse temperature (3 = 128f _1 , using 
the algorithm based on solution A. We have used the algorithm 
to study the quantum critical behavior of the one-dimensional 
Bose-Hubbard model with constant filling, using Renormal- 
ization Group flow equations. Studying the BKT transition is 
notoriously difficult because of the logarithmic finite-size cor- 
rections. The present algorithm has the big advantage of keep- 
ing the density constant, in contrast to the grand-canonical ap- 
proaches. The results of this study will be presented elsewhere 

We compare the efficiencies of solutions A and B. To this 
purpose, we simulate a one-dimensional lattice with 32 sites 
at an inverse temperature p = 32f ~' and a filling factor of one 
boson per site. The simulations consisted of 40 Markov chains 
that each ran 600 seconds after thermalization on a Pentium III 
processor. The same code was used, with only minor changes 
to go from algorithm A to B. We discuss the efficiency by 
looking at the standard deviations on the expectation value 
of V of Eq. d29i and on the average squared density. We 
calculated the squared density n 2 by averaging over all sites, 



(n 2 ) 



1 



X>*>- 



(33) 



The expectation value of the interaction term V was not cal- 
culated via the correlation function G(l), but by counting the 
number of interaction vertices in the configuration whenever 
the worm operator was diagonal 1 37] . When looking at the 
standard deviation on the squared density (Figure [3}, one can 
conclude solution A is the most efficient one. We found a sim- 
ilar picture when looking at the standard deviations on the ex- 
pectation value of Hq of Eq. ( 1291 . The errors on the standard 
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FIG. 2: Superfluid (p s ) and condensed fraction (p c ) for the one- 
dimensional homogeneous Bose-Hubbard model as a function of 
U //. The fractions have been calculated for a uniform lattice with 
N s = 128 at an inverse temperature (3 = 128? -1 , using Eqs. 13 1 1 and 
<32> . The condensed fraction was calculated from the correlation 
function (b\bi +r ), for which we have high statistics. 



deviations lie within ten percent. From the standard deviation 
on the expectation value of the interaction term V (Figure^), 
it follows that solution B is better in the superfluid phase. The 
same conclusion follows from the total energy. For the con- 
densed fraction the deviations are smallest for solution A for 
all values of U/t. Those on the superfluid fraction lie very 
close for solution A and B (see Figure[5}. We found that vary- 
ing the parameter (]) of Tables U and [n] does not change the 
efficiency in a significant way, as long as (j) is not too small. 
For further simulations we will always choose the parameter 
(j) as big as possible, under the constraint (j) < 9{j)d' ■ Figures 
|3j |4] and [5] show also standard deviations resulting from the 
directed loop SSE algorithm flrjl ITtL Il8il . One has to be very 
careful when comparing efficiencies of different algorithms. 
First, the SSE code works in the grand-canonical ensemble. 
In the SSE simulations, the chemical potential was changed 
in such a way that N remained constant. Second, the effi- 
ciency does not only depend on the algorithm, but also on the 
used data structures. In a SSE approach, the decomposition 
of the evolution operator corresponds to a perturbation expan- 
sion in all terms of the Hamiltonian, while the decomposition 
Eq. (0 perturbs only in the off-diagonal terms V. For the 
Bose-Hubbard model, where the contribution of the diagonal 
terms is large, the last approach is preferable. For all calcu- 
lated observables the standard deviations resulting from the 
SSE code were the largest. Figures[3]and|4]show the SSE de- 
viations increase rather rapidly with increasing U/t, whereas 
the deviations resulting from our method remain of the same 
order. We also calculated autocorrelation times for different 
observables. Here each bin ended after a constant number of 
measurements. We noticed that for solution B the autocor- 
relation times became very big for high values of U/t. For 
small U/t, the autocorrelation time for solution A is of the 
order of the number of Markov steps needed for 10 diagonal 



updates, and increases only slowly with increasing U/t. Of 
course it should be noted the algorithm based on solution A 
had to run much longer in order to get the same number of di- 
agonal measurements. For all measured observables we found 
similar autocorrelation times. 

We conclude that solution A, derived in the previous sec- 
tion, is more efficient than solution B, except in the super- 
fluid phase when looking at the interaction energy V. This 
can be understood by remarking that the time shifts of the 
worm operator are much larger for solution A. In the algo- 
rithm based on solution B, the time shifts are of the order of 
the mean imaginary time between two interaction vertices in 
the configurations. This explains why the standard deviations 
on the interaction energy are smallest for this solution in the 
superfluid phase. We also conclude that our algorithm is more 
efficient than the directed loop SSE algorithm when simulat- 
ing the one-dimensional Bose-Hubbard model. In the next 
section, we will apply the algorithm to a pairing model. In 
what follows, all calculations are performed using the algo- 
rithm based on solution B. 

0.01 f ■ ■ ' ' ' ■ ■ ' 3 

solution A 

solution B 

SSE 

0.001 - 



e 

1e-04 r 




1e-05 1 1 1 1 1 1 1 1 1 1 

1 23456789 10 

U/t 

FIG. 3: A comparison between the standard deviations on the 
squared density (see Eq. <33» resulting from the directed loop SSE 
algorithm and the algorithms based on solutions A and B. The homo- 
geneous Bose-Hubbard model is simulated for a lattice with 32 sites 
and 32 bosons at an inverse temperature (3 = 32t . The deviations 
resulted after a QMC calculation with 40 independent Markov chains 
that each ran 600 seconds on a Pentium III processor. 



IV. APPLICATION TO A PAIRING MODEL 

In the nuclear shell model, quantum Monte Carlo methods 
are valuable because they offer the possibility of doing calcu- 
lations in much larger model spaces than conventional diago- 
nalization techniques. Finite temperature shell model studies 
have been done with the aid of auxiliary-field QMC methods 
1231 . and ground-state properties of light nuclei have been cal- 
culated using variational and diffusion QMC techniques 1 38 ] . 
Furthermore, being able to calculate thermal properties of nu- 
clei makes it in principle possible to calculate nuclear level 
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FIG. 4: The standard deviation on the mean value of V (see Eq. 
1291 for the homogeneous Bose-Hubbard model as a function of U jt. 
Here solution A is the most efficient one in the Mott phase. In the 
superfluid phase, solution B becomes more efficient. 
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FIG. 5: The standard deviation on the mean value of the superfluid 
fraction p s . The deviations result from solutions A and B, and the 
directed loop SSE method. Each simulation consisted of 40 indepen- 
dent Markov chains that each ran 600 seconds. 



densities |39]. These densities are extremely important for 
making good theoretical estimates of nuclear reaction rates. 

The basic assumption in the shell model is the presence of a 
mean field in which the nucleons move. To improve on this, a 
residual interaction between the nucleons is introduced. Pair- 
ing between nucleons is the main short-range correlation in- 
duced by the residual interaction. Adding a simple pairing 
Hamiltonian to the mean-field Hamiltonian 



k. The operator n kt is the corresponding number-operator and 
e kt = e j, the single particle energy-eigenvalue. G t is the pair- 
ing strength for protons or neutrons. Proton-neutron pairing is 
not included, but this coupling contributes only in an impor- 
tant way for N = Z nuclei 1 42] . As a consequence, the prob- 
lem decouples for protons and neutrons. In the sequel only the 
neutron part of the model is considered, and the isospin index 
t is dropped to ease the notations. 

Based on algebraic techniques developed by Richardson 
lE^ll . the pairing model can be solved exactly for an arbitrary 
set of single particle levels at zero temperature lE4ll . In prac- 
tice it remains difficult to use these exact results to study the 
thermodynamics of the model, because the number of states 
needed in the ensemble increases very rapidly with increasing 
temperature. Thermodynamic al properties have been studied 
using auxiliary-field QMC, which is free of sign problems for 
the pairing Hamiltonian of Eq. J34l > when dealing with an even 
number of particles II24II . However, the present algorithm can 
consider nuclei with even and odd nucleon numbers. Note the 
auxiliary-field method scales as 0(N^) with N s the number 
considered single particle states, while a world-line algorithm 
scales linear with A^. 

When a nucleon occupies a single particle state k and its 
time-reversed state k is unoccupied, the nucleon is said to be 
'unaccompanied'. These states do not participate in the pair 
scattering by Hp. The mean-field plus pairing Hamiltonian 
can be rewritten as Eq. Q, 



H = H -V. 
H 



V = 



k L 

I fa 

H k^k' 



(35) 
(36) 

(37) 



The operators b\ — a\cij create a pair of nucleons in two time- 
reversed states and satisfy hard-core boson commutation rela- 
tions. In order to get the correct finite temperature properties, 
the possibility of changing the number of unaccompanied nu- 
cleons during the simulation should be incorporated. A path 
integral Monte Carlo method for the pairing Hamiltonian has 
been developed by Cerf and Martin |45J 14611 . but there the 
number of pairs remained fixed 1241 14711 . This problem can 
now be overcome elegantly by adding an extra pair breaking 
term 



V P eit = ^r£ E {b\a k ,a k n +H.c. 



(38) 



k k'^k" 



to the interaction part V of Eq. I37I . We define the worm 
operator as 



H, 



R P = £ ^k,n kt - £ -j Y t 4ft a W t a St a ^ W A = i E w * + J £ **** + f E E ( b l a k>ak» +H.c.) 



I 



mf - 



-p,n k 



t=p,n 



k.k' 



N ' 



' k^k' 



k k'^k" 



can account for this i4fl Elll . The operators at, create a 
particle in the single-particle eigenstate k of the mean-field 
Hamiltonian in the valence shell. The index t indicates pro- 
ton or neutron states and k denotes the time-reversed state of 



(39) 

with two extra parameters and g to be optimized. A term 
proportional to Vp el -t is included in the worm operator, in order 
to satisfy condition Eq. ( Il0i . This term will generate configu- 
rations with pair breaking interactions. However, it can occur 
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that too many of these interactions are generated, though we 
are only interested in generating configurations with a differ- 
ent number of unaccompanied particles, but without interac- 
tions of the type Eq. 08l >. This can be prevented by impos- 
ing the constraint that a configuration can contain at most two 
pair breaking interactions of this type. Observables are only 
updated if there are no pair-breaking interactions in the con- 
figuration. This means that a number of Markov steps are 
needed in order to reach a new allowed configuration with a 
different number of unaccompanied particles. When g of Eq. 
d39l is put equal to one, the percentage of diagonal configu- 
rations which contain no Vpert interactions (see Eq. J38» is 
about 15%. This is still efficient enough to sample the pairing 
Hamiltonian. There are a number of ways to increase the ef- 
ficiency. First of all one can change the parameter g, hereby 
influencing the appearance of pair breaking interactions. One 
can also restrict the number of times the worm tries to insert 
a Vpert interaction by allowing this only after a certain Markov 
time in which 'good' (i.e. without pair breaking interactions) 
configurations are sampled. One should also keep in mind that 
while a configuration contains pair breaking interactions, the 
worm itself is not necessarily of the pair breaking type. So a 
lot of Markov time is spend to change the configuration in a 
global way without removing the pair breaking interactions, 
leading to strong decorrelation. 

The main physical properties of nuclei in the Iron region 
are modeled by a schematic mean-field plus pairing Hamil- 
tonian. For the mean-field potential, we use a Woods-Saxon 
potential. Single particle energies are taken from Ref. II24I1 . A 
full fp + sdg valence space is chosen. These single-particle 
states and energies are shown in Table IIIII A pairing strength 
G = 16/56 MeV is used. Due to the size of the model space 
a strength smaller than the suggested value of 20 MeV per 
nucleon is used [24]. We have tested our code by compar- 
ing finite temperature results in a fp valence space with the 
ones obtained via an exact diagonalization technique H. We 
show results of calculations with the valence shell given in 
Table UlTl occupied by 10 and 11 valence neutrons. Figure [6] 
shows the expectation value of the neutron pairing-interaction 
operator (Hp) as a function of temperature. At low tempera- 
ture, the pairing energies are much lower for the even number 
of neutrons. This can be understood by remarking that for an 
odd number of neutrons there is always at least one unpaired 
nucleon. At temperatures higher than 1 MeV, the pairing ener- 
gies differ only slightly, because there is an increasing number 
of unpaired nucleons due to thermal excitation. This is also 
reflected in the specific heat (see Figure 0. A peak appears 
around 0.8 MeV due to the development of pair correlations. 

Because the worm operator conserves angular momentum, 
one can restrict the intermediate states to a specific value of 
the quantum numbers J and J z . This is not possible with the 
auxiliary-field QMC method. In our current implementation 
of the algorithm however, the occupation of each couple of 
time-reversed single-particle states (k,k) is exactly known at 
all times. Because the unaccompanied particle number opera- 
tor 

N" = Y^{n k -blb k ), (40) 



Single-particle energies (MeV) 


Orbital 


Protons 


Neutrons 


1/7/2 


-4.1205 


-10.4576 


2^3/2 


-2.0360 


-8.4804 


2P1/2 


1 'YX'XA 


- /.CO Iz 


1/5/2 


-1.2159 


-7.7025 


3*1/2 


4.7316 


-0.3861 


2^5/2 


5.6562 


0.2225 


2^3/2 


6.1324 


0.9907 


l 89/2 


6.6572 


0.5631 



TABLE III: Single particle eigenstates of a Woods-Saxon potential, 
taken from Ref. 12411 . The chosen valence space contains 42 states. 
The proton and neutron single particle energies (in MeV) are shown 
on the right. 



-1 




T[MeV] 

FIG. 6: Expectation value of the neutron part of the pairing- 
interaction operator as a function of temperature. The pairing 
strength G„ is equal to 16/56 MeV. We consider 10 and 11 neutrons 
in the fp + sdg valence space (see Tabel Hilt . At temperatures be- 
low 0.5 MeV the pairing energy is much lower for the even neutron 
number. 



commutes with the angular momentum projection operator J z 
(but not with 7 2 ), our current code allows restricting the sim- 
ulation to configurations with a fixed J z . Work on extending 
this technique to full J-projection is in progress. 

When the projection on J z was turned on, we included an 
extra global step in order to get a good convergence at the 
lowest temperatures. This extra global change allows for one 
or two unaccompanied nucleons (which block the state they 
occupy) to move to other states, and can occur whenever 
the worm is diagonal. First an unaccompanied nucleon at a 
blocked state I is chosen at random. A 'non-blocked' pair of 
states (k, k) is then chosen with probability 

P(k) = eJoMO+ttO-W/Ni, (41) 

with nk(t) the occupation number of state k at imaginary time 
f, and Ni a normalization factor. The subscript / indicates 
that the norm is determined for a configuration containing a 
blocked state I. The idea behind Eq. J41t is to get a probabil- 
ity distribution P(k) which is peaked around the Fermi level, 
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FIG. 7: The neutron specific heat C n as a function of temperature for 
10 and 1 1 neutrons in the full fp + sdg valence space (see Table lllH . 
Calculations were performed at a constant neutron pairing strength 
G„ = 16/56 MeV. 
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FIG. 8: J z projected internal energies as a function of temperature. 
The values of J z from to 7 are indicated on the left. A clear con- 
vergence to exact zero temperature results calculated via techniques 
explained in Ref. 1 44] can be seen. 



but other distributions can be chosen as well. The interchange 
of the occupations of the blocked pair of states (1,1) and the 
non-blocked pair (k, k) over the whole imaginary time interval 
P, is accepted with probability 

p=nrin(lA. (42) 

The acceptance factor for the case when the occupations of 
two pairs of non-blocked and blocked states are interchanged, 
can be constructed in a similar way. The extra step has a high 
acceptance rate, but is only necessary to enhance decorrela- 
tion at very low temperature when a J z -projection is included. 
At higher temperatures the unaccompanied nucleons move ef- 
ficiently from state to state via the last worm piece in Eq. d39l >. 
Figure [8] shows total energies after /..-projection at low tem- 
perature. Calculations were performed for ten neutrons mov- 
ing in the model space listed in Table [Hi] The neutron pair- 
ing strength is again G„ = 16/56 MeV. The figure also shows 
exact energy eigenvalues for J z — to J z — 7. These were 
calculated via a technique explained in Ref. 1 44] . The lowest 
J z = 1 , 2 and the lowest J z = 3 , 4 states are degenerate. For low 
enough values of T the finite temperature results clearly con- 
verge to the ground states within the considered ensembles. 

Note that we can compare with exact solutions because the 
pairing strength was taken constant for all levels. Our QMC 
method allows to solve pairing models with a single particle 
state dependent pairing strength G^y , for which no algebraic 
solutions are available. Taking in mind the method is appli- 
cable for even and odd nucleon systems and allows angular 
momentum symmetry projections, this could greatly extend 
the applicability of the pairing model. 

V. CONCLUSIONS AND OUTLOOK 

We have set up a quantum Monte Carlo method with a non- 
local loop updating scheme starting from a local worm op- 



erator in the path integral approach. This method allows to 
sample configurations with specific symmetries and, in partic- 
ular, to sample the canonical ensemble. It leads to a very ef- 
ficient sampling scheme with all moves accepted and without 
'bounces' or critical slowing down near second order phase 
transitions. We have proven detailed balance and tested er- 
godicity. Our method opens new perspectives for the study 
of quantum many-body systems where particle number and 
other symmetries play an important role. It can be applied 
to bosons, to fermions in absence of a sign problem and to 
non-frustrated spin systems at fixed magnetization. We have 
demonstrated this by simulating the Bose-Hubbard model and 
a nuclear pairing model. The equal-time one-body Green's 
function can be evaluated with high efficiency. When non- 
equal time observables are required, the current method can in 
principle still be combined with conventional non-local worm 
steps. There is still a lot of freedom in choosing the algorithm 
parameters, which can be used to optimize the algorithm. For 
the Bose-Hubbard model we compared the efficiency of our 
algorithm (with different parameter sets) with a directed loop 
SSE code. Though one should always be careful when com- 
paring different algorithms, we have strong indications that 
our method is very efficient. We have simulated a pairing 
model for even and odd particle numbers. Our finite temper- 
ature results clearly supplement algebraic methods and other 
QMC methods. Furthermore, a projection on angular momen- 
tum symmetries can be included. We have demonstrated this 
by showing J z -projected results. A work on full /-projection 
is in progress. 
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